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5.4 The Heat Equation and Convection-Diffusion 

The wave equation conserves energy. The heat equation u t = u xx dissipates energy. 
The starting conditions for the wave equation can be recovered by going backward in 
time. The starting conditions for the heat equation can never be recovered. Compare 
u t = cu x with ut = u xx , and look for pure exponential solutions u(x, t) = G(t) e lkx : 

Wave equation: G' = ickG G(t) = e ic ^ has \G\ = 1 (conserving energy) 
Heat equation: G' = —k 2 G G(t) = i has G < 1 (dissipating energy) 

Discontinuities are immediately smoothed out by the heat equation, since G is ex- 
ponentially small when k is large. This section solves u t = u xx first analytically and 
then by finite differences. The key to the analysis is the beautiful fundamental 
solution starting from a point source (delta function). We will show in equation (7) 
that this special solution is a bell-shaped curve: 

1 2 

u(x,t) = — -= e~ x I comes from the initial condition u(x,0) = 5(x) . (f) 



Notice that u t = cu x + du xx has convection and diffusion at the same time. The 
wave is smoothed out as it travels. This is a much simplified linear model of the 
nonlinear Navier-Stokes equations for fluid flow. The relative strength of convection 
by cu x and diffusion by du xx will be given below by the Peclet number. 

The Black-Scholes equation for option pricing in mathematical finance also has 
this form. So do the key equations of environmental and chemical engineering. 

For difference equations, explicit methods have stability conditions like At < 
|(Ax) 2 . This very short time step is more expensive than cAt < Ax. Implicit 
methods can avoid that stability condition by computing the space difference A 2 U 
at the new time level n + 1. This requires solving a linear system at each time step. 

We can already see two major differences between the heat equation and the wave 
equation (and also one conservation law that applies to both): 

1. Infinite signal speed. The initial condition at a single point immediately 
affects the solution at all points. The effect far away is not large, because of the 
very small exponential e~ x / 4 * in the fundamental solution. But it is not zero. 
(A wave produces no effect at all until the signal arrives, with speed c.) 

2. Dissipation of energy. The energy | J(u(x,t)) 2 dx is a decreasing function 
of t. For proof, multiply the heat equation u t = u xx by u. Integrate uu xx by 
parts with u(oo) = u{— oo) = to produce the integral of — (u x ) 2 : 

d Z* 00 ~\_ /*oo /»oo 

Energy decay — / —u 2 dx= / uu xx dx = — / (u x ) 2 dx < . (2) 

dt J -oo 2 J -oo J -oo 
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3. Conservation of heat (analogous to conservation of mass): 



Heat is conserved 



d rx 
dt 



u(x, t) dx 



u x (x,t) 



(3) 



J x=— oo 



Analytic Solution of the Heat Equation 
Start with separation of variables to find solutions to the heat equation: 

Ql fill 

Assume u(x,t) = G(t)E(x). Then u t = u xx gives G E = GE" and — = — — . 

G E 

(4) 

The ratio G'/G depends only on t. The ratio E" j 'E depends only on x. Since 
equation (4) says they are equal, they must be constant. This produces a useful 
family of solutions to u t = u xx . 



E" G' ■> 
—— = — is solved by E(x) = e lkx and G(t) = e~ k " 1 

E G 



Two x- derivatives produce the same — k 2 as one t-derivative. We are led to exponential 
solutions of e tkx e~ k t and to their linear combinations (integrals over different k): 

1 f°° 2 

General solution u(x,t) = — / u (k)e lkx e~ k t dx. (5) 

2tt J _ 



-oo 



At t — 0, formula (5) recovers the initial condition u(x, 0) because it inverts the 
Fourier transform uq (Section 4.4.) So we have the analytical solution to the heat 
equation — not necessarily in an easily computable form ! This form usually requires 
two integrals, one to find the transform Uo(k) of u(x,0), and the other to find the 
inverse transform of Uo(k)e~ h 1 in (5). 

Example 1 Suppose the initial function is a bell-shaped Gaussian u(x,0) = e~ x2 l 2a . 
Then the solution remains a Gaussian. The number a that measures the width of the 
bell increases to a + 2t at time t, as heat spreads out. This is one of the few integrals 

2 

involving e~ x that we can do exactly. Actually, we don't have to do the integral. 

That function e _a;2//2cr is the impulse response (fundamental solution) at time t = 
to a delta function S(x) that occurred earlier at t = So the answer we want (at 

time t) is the result of starting from that 5(x) and going forward a total time \o + 1: 



Widening Gaussian u(x, t) = e -^ 2 /(^ + ^) . ( 6 ) 



This has the right start at t = and it satisfies the heat equation. 
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The Fundamental Solution 

For a delta function u(x, 0) = 5(x) at t = 0, the Fourier transform is Uo(k) = 1. Then 
the inverse transform in (5) produces u(x,t) = f e lkx e~ k 1 dk One computation of 
this u uses a neat integration by parts for du/dx. It has three — l's, from the integral 
of ke~ k t and the derivative of ie lkx and integration by parts itself: 



This linear equation du/dx = —xu/2t is solved by u = ce~ x2 ^ u . The constant 
c = 1 / V 47rt is determined by the requirement J u(x,t)dx = 1. (This conserves 
the heat J it(x, 0) dx = J 5(x) dx = 1 that we started with. It is the area under a 
bell-shaped curve.) The solution (1) for diffusion from a point source is confirmed: 

Fundamental solution from , s 1 —rAlAi 

( n\ xr \ u(x,t) = —==e x I* 1 . (8) 

u(x,0) = S(x) v ^rt 



In two dimensions, we can separate x from y and solve u t = u xx + u yy : 

u(x,y,t) = [ ^= ) e -*-/4* e -ir/«. ( 9 ) 



2 

Fundamental solution from , +\ — ( ^ \ -% 2 /4* /4* 

w(a;,?/,0) = S(x)S(y) 

With patience you can verify that u(x,t) and u(x,y,t) do solve the ID and 2D heat 

equations (Problem ). The zero initial conditions away from the origin are 

correct as t — > 0, because e~ c /* goes to zero much faster than \/\ft blows up. And 
since the total heat remains at J u dx = 1 or J J udx dy = 1, we have a valid solution. 

If the source is at another point x — s, then the response just shifts by s. The 
exponent becomes — (x— s) 2 /At instead of — x 2 /At. If the initial u(x, 0) is a combination 
of delta functions, then by linearity the solution is the same combination of responses. 
But every u(x, 0) is an integral J 5(x — s) u(s, 0) ds of point sources ! So the solution to 
ut = u xx is an integral of the responses to 8(x — s). Those responses are fundamental 
solutions starting from all points x = s: 



Solution from any u(x, 0) u(x,t) = / u(s, 0) e~( x ~ s ) / At ds . (10) 




Now the formula is reduced to one infinite integral — but still not simple. And for a 
problem with boundary conditions at x = and x = 1 (the temperature on a finite 
interval, much more realistic), we have to think again. Similarly for an equation 
Ut = (c(x)u x ) x with variable conductivity or diffusivity. That thinking probably 
leads us to finite differences. 

I see the solution u(x, t) in (10) as the convolution of the initial function u(x, 0) 
with the fundamental solution. Three important properties are immediate: 
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1. // u(x, 0) > for all x then u(x,t) > for all x and t. Nothing in 
formula (10) will be negative. 

2. The solution is infinitely smooth. The Fourier transform uo(k) in (5) is 
multiplied by e~ k 1 . In (10), we can take all the x and t derivatives we want. 

3. The scaling matches x 2 with t. A diffusion constant d in the equation 
u t = du xx will lead to the same solution with t replaced by dt, when we write 
the equation as du/d{dt) = d 2 u/dx 2 . The fundamental solution has e~ x ' 4 * 
and its Fourier transform has e~ dk 1 . 



Example 2 Suppose the initial temperature is a step function u(x, 0) = 0. Then for 
negative x and u(x, 0) = 1 for positive x. The discontinuity is smoothed out immediately, 
as heat flows to the left. The integral in formula (10) is zero up to the jump: 

u(x, t) = -== / e< x ~ s ) l At ds . (11) 

No luck with this integral ! We can find the area under a complete bell-shaped curve 
(or half the curve) but there is no elementary formula for the area under a piece of the 

2 

curve. No elementary function has the derivative e~ x . That is unfortunate, since those 
integrals give cumulative probabilities and statisticians need them all the time. So they 
have been normalized into the error function and tabulated to high accuracy: 

2 f x _ 2 

Error function erf (x) = —= I e s ds . (12) 

V ^ Jo 

The integral from — x to is also erf(x). The normalization by 2/y/n gives erf(oo) = 1. 

We can produce this error function from the heat equation integral (11) by setting 
S = (s — x)/\fM. Then s = changes to S = —x/\^4t as the lower limit on the integral, 
and dS = ds/y/At. Split into an integral from to oo, and from —x/\/At to 0: 

u(x,t) = ^E = [°° e~ s2 dS = -(l + erf ( ^=Y) . (13) 

Good idea to check that this gives u = \ at x = (where the error function is zero). 
This is the only temperature we know exactly, by symmetry between left and right. 



Explicit Finite Differences 

The simplest finite differences are forward for du/dt and centered for d 2 u/dx 2 : 

Explicit method ^ = U ^\~ ^ = = U ^ n . (14) 

P At (Ax) 2 At (Ax) 2 v ; 
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Each new value L^, n +i is given explicitly by Uj >n + R(Uj + i in — 2Uj^ n + Uj^i). The 
crucial ratio for the heat equation ut = u xx is now R = At /(Ax) 2 . 

We substitute U^ n = G n e ikjAx to find the growth factor G = G(k, At, Ax): 

One-step Q=l + R(e ikAx -2 + e ~ lkAx ) = 1 + 2i*(cos fcA^ - 1) . (15) 

growth factor 

G is real, just as the exact one-step factor e _fc2Ai is real. Stability requires \G\ < 1. 
Again the most dangerous case is when the cosine equals —1 at kAx = n: 

At 1 

Stability condition \G\ = 1 — AR\ < 1 which requires R = < — . (16) 

(Aa;) 2 2 

In many cases we accept that small time step At and use this simple method. The 
accuracy from forward A t and centered A 2 X is \U — u\ = 0(At + (Ax) 2 ). Those two 
error terms are comparable when R is fixed. 

We could improve this one-step method to a multistep method. The "method 
of lines" calls an ODE solver for the system of differential equations (continuous in 
time, discrete in space). There is one equation for every meshpoint x = jh: 

Method of Lines = — - — — - 1 = 3+ — \ n J — - . (17) 

dt (Ax) 2 dt (Ax) 2 y 1 

This is a stiff system, because its matrix —K (second difference matrix) has a large 
condition number: A max (i^)/A min (i^) ~ N 2 . We could choose a stiff solver like odel5s 
in MATLAB. 



Implicit Finite Differences 

A fully implicit method for u t = u xx computes A 2 U at the new time (n + l)Ai: 

t i * • i A. t U n A^Un+i Uj n+ i—Uj n Uj + \ n+ \ — 2Uj n+ \ + Uj-i n +l , 1 oN 

Implicit = — — — — = — — — . 18 

H At (Ax) 2 At (Ax) 2 v ' 

The accuracy is still first-order in time and second-order in space. But stability no 
longer depends on the ratio R = At /(Ax) 2 . We have unconditional stability, with a 
growth factor < G < 1 for all k. Substituting Uj >n = G n e l ^ kAx into (18) and then 
canceling those terms from both sides leaves an extra G on the right side: 

G = 1 + RG(e lkAx - 2 + e ~ lkAx ) leads to G= — — - . (19) 

v ; 1 + 2R(1 - cos kAx) v ' 

The denominator is at least 1, which ensures that < G < 1. The time step is 
controlled by accuracy, because stability is no longer a problem. 
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There is a simple way to improve to second-order accuracy. Center everything at 
step n + ~. Average an explicit A 2 x U n with an implicit A 2 U n+ i. This produces the 
famous Crank- Nicols on method (like the trapezoidal rule): 

Crank-Nicolson Uj ' n+ \~ Uj ' n = —^(A 2 x U j>n + A 2 x U j>n+1 ) . (20) 



Now the growth factor G, by substituting Uj >n = G n e^ kAx into (20), solves 

G — 1 (jr "I - 1 , , , s 

— — = — — — 2 cos k Ax - 2 . 21 

At 2(Ax) 2 v ; v ' 

Separate out the part involving G, write R for Ai/(Ax) 2 , and cancel the 2's: 

Unconditional stability G = 7 — f has \G\ < 1 . (22) 

1 - i?(cos A;Ax - 1) 1 1 - v ; 

The numerator is smaller than the denominator, since cos Arc < 1. We do notice 
that cos A; Ax = 1 whenever kAx is a multiple of 2ir. Then G = 1 at those frequencies, 
so Crank-Nicolson does not give the strict decay of the fully implicit method. We 
could weight the implicit A 2 JJ n+ i by a > | and the explicit A 2 x U n by 1 — a < |, to 
give a whole range of unconditionally stable methods (Problem ). 



Numerical example 



Finite Intervals with Boundary Conditions 

We introduced the heat equation on the whole line — 00 < x < 00. But a physical 
problem will be on a finite interval like < x < 1. We are back to Fourier series 
(not Fourier integrals) for the solution u(x, t). And second differences bring back the 
great matrices K,T,B,C that depend on the boundary conditions: 

Absorbing boundary at x = 0: The temperature is held at u(0, t) = 0. 
Insulated boundary: No heat flows through the left boundary if 1^(0, t) = 0. 

If both boundaries are held at zero temperature, the solution will approach u(x, t) = 
everywhere as t increases. If both boundaries are insulated as in a freezer, the solution 
will approach u(x,t) = constant. No heat can escape, and it is evenly distributed as 
t — > 00. This case still has the conservation law u(x, t) dx = constant. 

Example 3 (Fourier series solution) We know that e tkx is multiplied by e _fc2 * to give 
a solution of the heat equation. Then u = e~ h2t sinkx is another solution (combining 
+k with —k). With zero boundary conditions w(0,t) = u(l,t) = 0, the only allowed 
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frequencies are k = nn (then sinnnx = at both ends x = and x = 1). The complete 
solution is a combination of these exponential solutions with k = nn: 

oo 

Complete solution u(x, t) = b n e~ n n t sin mix . (23) 

n=l 

The Fourier sine coefficients b n are chosen to match u(x, 0) = ^6 n sinn7ra; at t = 0. 

You can expect cosines to appear for insulated boundaries, where the slope (not 
the temperature) is zero. This gives exact solutions to compare with finite difference 
solutions. For finite differences, absorbing boundary conditions produce the matrix 
K (not B or C). The choice between explicit and implicit decides whether we have 
second differences —KU at time level n or level n + 1: 

Explicit method U n+1 - U n = —RKU n (24) 

Fully implicit U n+1 -U n = -RKU n+1 (25) 

Crank-Nicolson U n+l -U n = -\RK(U n + U n+1 ) . (26) 

The explicit stability condition is again R < | (Problem ). Both implicit meth- 
ods are unconditionally stable (in theory). The reality test is to try them in practice. 

An insulated boundary at x = changes K to T. Two insulated boundaries 
produce B. Periodic conditions will produce C. The fact that B and C are singular 
no longer stops the computations. In the fully implicit method (/ + RB)U n+ \ = U n , 
the extra identity matrix makes I + RB invertible. 

The two-dimensional heat equation describes the temperature distribution in 
a plate. For a square plate with absorbing boundary conditions, the difference matrix 
K changes to K2T). The bandwidth jumps from 1 (triangular matrix) to N (when 
meshpoints are ordered a row at a time). Each time step of the implicit method 
now requires a serious computation. So implicit methods pay an increased price for 
stability, to avoid the explicit restriction At < |(Ax) 2 + |(Ay) 2 . 



Convection-Diffusion 

Put a chemical into flowing water. It diffuses while it is carried along by the flow. A 
diffusion term du xx appears together with a convection term cu x . This is the simplest 
model for one of the most important differential equations in engineering: 

du du d^u 
Convection-diffusion equation — = c — — h d 7—^ . (27) 

at ox ox 1 

On the whole line — 00 < x < 00, the flow and the diffusion don't interact. If the 
velocity is c, convection just carries along the diffusing solution to h t = dh x 



''XX ■ 



Diffusing traveling wave 



u(x, t) — h(x + ct, t) . 



(28) 
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Substituting into equation (27) confirms that this is the solution (correct at t = 0): 

. du dh dh dh , d 2 h du , d 2 u , , 

Cham rule — = c— + — = c— + d -— = c— + d — - . (29) 

dt ox dt ox ox 1 ox ox 1 

Exponentials also show this separation of convection e lkct from diffusion e~ dk2t : 
Starting from e ikx u(x, t) = e~ dkH e ik ( x + c *> . (30) 

Convection-diffusion is a terrific model problem, and the constants c and d clearly 
have different units. We take this small step into dimensional analysis: 

„ „ distance (distance) 2 

Convection coefficient c: ; Diffusion coefficient a: ; (31) 

time time 

Suppose L is a typical length scale in the problem. The Peclet number Pe = cL/d 
is dimensionless. It measures the relative importance of convection and diffusion. This 
Peclet number for the linear equation (27) corresponds to the Reynolds number for 
the nonlinear Navier-Stokes equations (Section ). 

In the difference equation, the ratios r = cAt/Ax and 2R = 2dAt/(Ax) 2 are also 
dimensionless. That is why the stability conditions r < 1 and 2R < 1 were natural for 
the wave and heat equations. The new problem combines convection and diffusion, 
and the cell Peclet number P uses Ax/2 as the length scale in place of L: 

Cell Peclet Number P = — = . (32) 

2R 2d K J 

We still don't have agreement on the best finite difference approximation! Here 
are three natural candidates (you may have an opinion after you try them): 

1. Forward in time, centered convection, centered diffusion 

2. Forward in time, upwind convection, centered diffusion 

3. Explicit convection (centered or upwind), with implicit diffusion. 

Each method will show the effects of r and R and P (we can replace r/2 by RP): 

1. Centered explicit ^' w+1 ~ Uj ' n = c Uj+1 ' n ~ + d . (33) 

A* 2Ax (Ax) 2 v ; 

Every new value Uj >n +i is a combination of three known values at time n: 

U j>n+1 = (1 - 2R)U hn + (R+ RP)U j+1>n + (R— RP)U^ 1>n . (34) 

Those three coefficients add to 1, and U = constant certainly solves equation (33). // 
all three coefficients are positive, the method is surely stable. More than 
that, oscillations cannot appear. Positivity of the middle coefficient requires R < |, 
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as usual for diffusion. Positivity of the other coefficients requires |P| < 1. Of course 
P will be small when Ax is small (so we have convergence as Ax — > 0). In avoiding 
oscillations, the actual cell size Ax is crucial to the quality of U. 

Figure 5.12 was created by Strikwerda [59] and Persson to show the oscillations for 
P > 1 and the smooth approximations for P < 1. Notice how the initial hat function 

is smoothed and spread and shrunk by diffusion. Problem finds the exact 

solution, which is moved along by convection. Strictly speaking, even the oscillations 
might pass the stability test \G\ < 1 (Problem ). But they are unacceptable. 



Figure 5.12: Convection-diffusion with and without numerical oscillations: R = 
, r = and . 



2tt • J j.- Uj n +\ Uj n Uj+\ n Uj n A Uj n /nr\ 

. Upwind convection — — J -^ = c— — — V d , A ' . (35) 

At Ax (Ax) 2 v ; 

The accuracy in space has dropped to first order. But the oscillations are eliminated 
whenever r + 2R < 1. That condition ensures three positive coefficients when (35) is 
solved for the new value Uj, n+ i. 

U jin+1 = (RP + R)U j+1 , n + (1 - RP - 2R)U j>n + RUj-^n ■ (36) 

Arguments are still going, comparing the centered method 1 and the upwind method 2. 
The difference between the two convection terms, upwind minus centered, is ac- 
tually a diffusion term hidden in (35) ! 

Extra diffusion U '*\ = U -i - ^ = ^ = (**\ U '» ~-2^± ^ri . (37) 



Ax 2Ax V 2 / ( Ax ) 2 

So the upwind method has this extra numerical diffusion or "artificial viscosity" 
to kill oscillations. It is a non-physical damping. If the upwind approximation were 
included in Figure 5.12, it would be distinctly below the exact solution. Nobody is 
perfect. 

3. Implicit diffusion U ^+\~ U ^ = c U i+i,n -U j>n A x*>Wi (3g) 

At Ax (Ax) 2 v ; 

MORE TO DO 

Problem Set 5.4 

1 Solve the heat equation starting from a combination u(x, 0) = 5(x + 1) — 25(x) + 
S(x — 1) of three delta functions. What is the total heat J u(x, t) dx at time t ? 
Draw a graph of u(x, 1) by hand or by MATLAB. 

2 Integrating the answer to Problem 1 gives another solution to the heat equation: 

Show that w(x,t) = / u(X,t)dX solves w t = w xx . 

Jo 

Graph the initial function w(x, 0) and sketch the solution w(x, 1). 
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3 Integrating once more solves the heat equation h t = h xx starting from h(x, 0) = 
/ w(X, 0) dX = hat function. Draw the graph of h(x, 0). Figure 5.12 shows the 
graph of h(x,t), shifted along by convection to h(x + ct,t). 

4 In convection-diffusion, compare the condition R < |,P < 1 (for positive coef- 
ficients in the centered method) with r + 2R < 1 (for the upwind method). For 
which c and d is the upwind condition less restrictive, in avoiding oscillations? 

5 The eigenvalues of the n by n second difference matrix K are = 2 — 2 cos 
The eigenvectors y k in Section 1.5 are discrete samples of smknx. Write the 
general solutions to the fully explicit and fully implicit equations (14) and (18) 
after N steps, as combinations of those discrete sines y k times powers of 



Another exact integral involving e x ' is 



_ „2 

x e 



x2 l u dx 



-2te 



x 2 /U 



It. 



From (17), show that the temperature is u = \ft at the center point x = 
starting from a ramp u(x, 0) = max(0, x). 

A ramp is the integral of a step function. So the solution of u t = u xx starting 
from a ramp (Problem 6) is the integral of the solution starting from a step 
function (Example 2 in the text). Then \rt must be the total amount of heat 
that has crossed from x > to x < in Example 2 by time t. Explain each of 
those three sentences. 



